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I. INTRODUCTION 



Nuclear structure models based on energy density functionals (EDFs) have successfully 
been used over the whole nuclide chart, from relatively light systems to superheavy nuclei, 
and from the valley of /5-stability to the particle drip-lines . In lowest order - the 

mean-field approximation, an EDF is constructed as a functional of one-body nucleon den- 
sity matrices that correspond to a single product state - Slater determinant of single-particle 
or single-quasiparticle states. This framework can thus also be referred to as single reference 
(SR) EDF. The static nuclear mean-field is characterized by symmetry breaking - transla- 
tional, rotational, particle number. Even though symmetry breaking incorporates important 
static correlations, i.e. deformations and pairing in the SR EDF, this framework can only 
describe ground-state properties: binding energies, charge radii, etc. Excitation spectra 
and electromagnetic transition probabilities can only be calculated by including correlations 
beyond the static mean-field through restoration of broken symmetries and configuration 
mixing of symmetry-breaking product states. The most effective approach to configuration 
mixing calculations is the generator coordinate method (GCM) {4], with multipole moments 
used as collective coordinates that generate the symmetry-breaking product wave functions. 
In such a multi-reference (MR) EDF approach, families of static mean-field configurations 
are mixed to restore symmetries and take into account fluctuations of collective variables. 
The corresponding EDFs are functionals of transition densities built from pairs of symmetry- 
breaking product states. 

In the first two parts of this work we have extended the framework of relativistic 

energy density functionals to include correlations related to the restoration of broken sym- 
metries and to fluctuations of collective variables. A model has been developed that uses the 
GCM to perform configuration mixing of angular-momentum and particle-number projected 
relativistic wave functions. The geometry is restricted to axially symmetric shapes, and the 
intrinsic wave functions are generated from the solutions of the relativistic mean-field + 
Lipkin-Nogami BCS equations, with a constraint on the mass quadrupole moment. The 
model employs a relativistic point-coupling (contact) nucleon-nucleon effective interaction 
in the particle-hole channel, and a density-independent ^-interaction in the particle-particle 
channel. This approach enables a quantitative description of the evolution of shell-structure, 
deformation and shape coexistence phenomena in nuclei with soft potential energy surfaces. 
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In the first application [3] , the GCM based on relativistic EDFs was employed in a study 
of shape transitions in Nd isotopes. It has been shown that the microscopic framework based 
on universal EDFs, adjusted to nuclear ground-state properties, and extended to take into 
account correlations related to symmetry restoration and fluctuations of collective variables, 
describes not only general features of shape transitions, but also the unique behavior of the 
excitation spectra and transition rates at the X(5) critical point of quantum shape phase 
transition. However, an exact diagonalization of the X(5) Hamiltonian carried out in Ref. j^, 
has shown that many properties of the solution are dominated by /5 — 7 coupling induced by 
the kinetic energy operator. The importance of the explicit treatment of the triaxial degree 
of freedom, i.e. inclusion of /? — 7 coupling, was also emphasized in two recent studies of 



shape transitions in the Nd isotopic chain 



10|, that used the self-consistent Hartree-Fock- 



Bogoliubov model, based on the finite-range and density-dependent Gogny interaction, to 
generate potential energy surfaces in the /? — 7 plane. 

While GCM configuration mixing of axially symmetric mean-field states has been imple- 
mented by several groups and routinely used in nuclear structure studies, the apphcation of 
this method to triaxial shapes is a much more difficult problem. Only very recently a model 
has been introduced based on the mean-field states generated by triaxial quadrupole 
constraints that are projected on particle number and angular momentum and mixed by 
the generator coordinate method. This method is equivalent to a seven-dimensional GCM 
calculation, mixing all five degrees of freedom of the quadrupole operator and the gauge an- 
gles for protons and neutrons. However, the numerical implementation of the model is very 
complex, and applications to medium-heavy and heavy nuclei are still computationally too 
demanding and time consuming. In addition, the use of general EDFs, i.e. with an arbitrary 
dependence on nucleon densities, in GCM type calculations, often leads to discontinuities 
or even divergences of the energy kernels as a function of deformation [121, 13] • Only for 
certain types of density dependence a regularization method can be implemented [l^, which 
corrects energy kernels and removes the discontinuities and divergences. 

In an alternative approach to five-dimensional quadrupole dynamics that includes rota- 
tional symmetry restoration and takes into account triaxial quadrupole fiuctuations, a collec- 
tive Bohr Hamiltonian is constructed, with deformation-dependent parameters determined 
from microscopic self-consistent mean-field calculations [isIJigI. The collective Hamiltonian 
can be derived in the Gaussian overlap approximation (GOA) [4] to the full five- dimensional 
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GCM. With the assumption that the GCM overlap kernels can be approximated by Gaus- 
sian functions, the local expansion of the kernels up to second order in the non-locality 
transforms the GCM Hill- Wheeler equation into a second-order differential equation - the 
Schrodinger equation for the collective Hamiltonian. The kinetic part of this Hamiltonian 
contains an inertia tensor 1^], and the potential energy is determined by the diagonal ele- 
ments of the Hamiltonian kernel, and also includes zero-point energy (ZPE) corrections 18|. 
The adiabatic time-dependent Hartree-Fock (ATDHF) theory 19| provides an alternative 
way to derive a classical collective Hamiltonian and, after re-quantization, a Bohr Hamil- 
tonian of the same structure is obtained, but with different microscopic expressions for the 



mertia parameters 
the collective Bohr Hamiltonian 
Yoccoz masses 



201 ] ■ There is a long-standing debate in the literature about masses in 



23|, i.e. whether the GCM-GOA expressions (the so-called 



2l|) or the ATDHF expressions (the so-called Thouless-Valatin masses 22l |) 



should be used. The Thouless-Valatin masses have the advantage that they also include the 
time-odd components of the microscopic wave functions and, in this sense, the full dynamics 
of a nuclear system. In the GCM approach these components can only be included if, in 
addition to the coordinates g^, the corresponding canonically conjugate momenta pi are also 
taken into account, but this is obviously a very complicated task. In many applications 
a further simplification is thus introduced in terms of cranking formulas [18|, l2J], i.e. the 
perturbative limit for the Thouless-Valatin masses, and the corresponding expressions for 
ZPE corrections. This approximation was applied in recent studies using models based both 
on the Gogny interaction 25|], and Skyrme energy density functionals 261]. 

In this work we develop a new implementation for the solution of a five- dimensional col- 
lective Hamiltonian that describes quadrupole vibrational and rotational degrees of freedom, 
with parameters determined in the framework of relativistic EDF. An initial study along 
this line which, however, did not include ZPE corrections, was reported in Ref. j27| . 

The theoretical framework is described in Sec.[Tll the method of solution of the eigenvalue 
problem of the general collective Hamiltonian, and the calculation of the mass parameters, 
moments of inertia, and ZPE corrections. In Sec. Illll the model is tested in the calculation of 
collective excitation spectra of the chain of even-even Gd isotopes, and results are compared 
with available data. Sec. [IV]presents a summary and an outlook for future studies. Technical 
details about the solution of the Dirac equation in triaxial geometry, the calculation of 
moments of inertia, ZPE corrections, and numerical tests are included in Appendix A-D. 
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II. THEORETICAL FRAMEWORK 



A. Collective Hamiltonian in five dimensions 

Nuclear excitations determined by quadrupole vibrational and rotational degrees of 
freedom can be treated simultaneously by considering five quadrupole collective coordi- 
nates a^, n = —2,-1,..., 2 that describe the surface of a deformed nucleus: R = 
Rq [1 + ^fj.^2fi] ■ To separate rotational and vibrational motion, these coordinates are 
usually parameterized in terms of two deformation parameters f3 and 7, and three Euler 
angles (0, 6, ip) = Q which define the orientation of the intrinsic principal axes in the 
laboratory frame 

a, = D%{n)[3cos^ + ^ [D%in) + Dl_,in)] /3sin7 , (1) 

where D^j, is the Wigner function [28]. The three terms of the classical collective Hamilto- 
nian, expressed in terms of the intrinsic variables /3, 7 and Euler angles 

Hcon = Xih{f3, 7) + ^rot(/3, 7, ^) + Vcoii(/5, 7) , (2) 
denote the contributions from the vibrational kinetic energy: 

Tvib = ^B^pP' + ISBp^h + \l3^B,,^^ , (3) 



the rotational kinetic energy: 

1 ^ 

^-t = -5ZXfc^2^ (4) 



fc=i 

and the collective potential energy Vcoii(/5, 7)- The mass parameters Bp^^, Bp^, B^^, and the 
moments of inertia depend on the quadrupole deformation variables (3 and 7. 



The Hamiltonian Eq. ([2]) is quantized according to the general Pauli prescription |29| : 
for the classical kinetic energy 

T = ]^^Bij{q)qiqj , (5) 
and the corresponding quantized form reads: 



kin 



2 Vdet5 ^ dqi dqj 
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The kinetic energy tensor in Eq. ([2]) takes the block diagonal form: 



B 



with the vibrational part of the tensor 



-Bvib 



B 



vib 







B 



(7) 



rot 



Bi3f3 /3B/3^ 

f3Bf^<y /3'^By'y 



(8) 



In general the rotational part is a complicated function of the Euler angles but, using the 
quasi-coordinates related to the components of the angular momentum in the body-fixed 
frame, it takes a simple diagonal form 



8ikXk, — 1, 2, 3 



(9) 



with the moments of inertia expressed as 



Ik = 4:Bkf3^ sin^(7 - 2kn/3) 



(10) 



This particular functional form is motivated by the fact that all three moments of inertia 
vanish for the spherical configuration {f3 = 0) and, additionally, Xz and Xy vanish for axially 
symmetric prolate (7 = 0*^) and oblate (7 = 60°) configurations, respectively. The resulting 
determinant reads 

detB = det^^ib ■ det^^ot = 4wr/5^ sin^ 87 , (11) 

where w = BppB^^ — B'^^^ and r = B1B2B3. The quantized collective Hamiltonian can hence 
be written in the form: 



H — Tvib + Trot + Vc 



coll 



(12) 



with 



- vib 



+ 



2./wr {(3"^ 
1 



d rr 



d d fr 



8(3 \ w 



dp dpy w 



P sin 87 



d rr 



sin 3'jB/s^ 



d 1 d [r 
+ 



sin 875^/3 



d 



and 



d'y\w~' " ' ^ 8(3 ' /3 d'y \ w ~ " ' ^7 
P 



- rot 



1 ^ 

E 



(18) 



(14) 



k=l 
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where Jk denotes the components of the angular momentum in the body-fixed frame of a 
nucleus. Vcow is the collective potential. The Hamiltonian describes quadrupole vibrations, 
rotations, and the coupling of these collective modes. The determinant Eq. (fTT!) determines 
the volume element in the collective space: 

-2^ 



dr. 



coll 



dildroy/wr 



dj3(3^ / (i7|sin37| / d^l^/wr . 
Jo 



(15) 



and the quantized Hamiltonian Eq. (|T2|) is hermitian with respect to the collective measure 
Eq. (USD. 

The methods used to solve the eigenvalue problem of the general collective Hamiltonian 
Eq. (fT2l) can be divided into two classes. The first is based on a direct numerical solution 
of a system of partial differential equations using finite-difference methods [30|, l3l|, l32] • The 
second approach uses an expansion of eigenf unctions in terms of a complete set of basis 
functions, that depend on the deformation variables P and 7, and the Euler angles 0, 9 



and ip [23 



34 



35 



36[. The eigenvalue problem reduces to a simple matrix diagonalization. 



and the main task is the construction of an appropriate basis for each value of the angular 
momentum quantum number. 



In this work we employ t 
method described in Refs. 



le second approach and construct basis states according to the 



25 



36 



37 



38 



39j. For each value of the angular momentum /, 



one chooses a complete set of square integrable functions 



cosm7 
sin rwy 



ML 



(16) 



The projections M and L are determined by the angular momentum: M, L = — /, ...,/. In 
principle, the parameter n can take any nonnegative integer value, but in actual calculations 
a certain cut-off value rimax has to be imposed. The allowed values of m are: m = n,n — 
2, . . . , or 1. The choice of the function e"'^^^^/^ ensures that the basis states generate wave 
functions that vanish at large deformations {(3 00). The basis parameter /i has to be 
adjusted for each nucleus individually, so that it minimizes the ground state energy of the 
nucleus. However, if the cut-off value rimax is large enough, a stable ground-state solution 
can be found for a broad range of values of the parameter /i . 

The basis states have to fulfill certain symmetry conditions that originate from the fact 
that the choice of the body-fixed frame is not unique. For a given quadrupole tensor 
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in the laboratory frame, there are 24 possible orientations of the body-fixed right-hand 
coordinate system, corresponding to different values of the variables jS, 7, and Q. The basis 
states in the body-fixed frame must be invariant with respect to the transformations that 
connect various choices of the body-fixed frame, and which form a finite group isomorphic 

3C 



to the octahedral point group O 



40| (group of proper rotations which take a cube or 



octahedron into itself). This symmetry condition is fulfilled by linear combinations of the 
states ffT6l) 



KeAI 



(17) 



invariant under the transformations of the octahedral group. The angular part corresponds 
to linear combinations of the Wigner functions 



21 + 1 



and the summation in Eq. (1171) is over the allowed set of the K values: 



(18) 



A/ 



0,2, 



/ for / mod 2 = 



2,4 / - 1 for I mod 2 



(19) 



In the next step linearly independent functions have to be selected from the overcomplete 
basis set Eq. f|T7j) . In addition, some of the basis states have to be discarded in order 
to enforce the correct behavior of solutions on the 7 = nir/S axes 30|]. A simple and 
elegant solution of both problems is provided by group theoretical methods 4l||. Finally, 
the basis states Eq. fll7l) are not orthogonal. Although the Hamiltonian could also be 
diagonalized directly in a non-orthogonal basis [35|, we choose to orthogonalize the basis 



42| 



states by applying the Cholesky-Banachiewicz method 

The diagonalization of the collective Hamiltonian yields the energy spectrum and the 
corresponding eigenfunctions 



(20) 



KeAi 



Using the collective wave functions Eq. (!20l) . various observables can be calculated and 
compared with experimental results. For instance, the quadrupole E2 reduced transition 
probability: 



B{E2; al a' I' 



1 



2/ + 1 
8 



\{a'I'\\M{m\<^I)\' 



(21) 



and the spectroscopic quadrupole moment of the state \al): 

Qspecai = -^^==CjUc^I\\MiE2)\\aI) , (22) 
where jCi {E2) denotes the electric quadrupole operator. Detailed expressions for the reduced 



matrix element {a' I'\\Ai{E2)\\aI) can be found in Ref. 30|. 

The shape of a nucleus can be characterized in a qualitative way by average values of the 
invariants (3^ cos 87, as well as their combinations. For example, the average value of the 
invariant in the state \al): 

{f3')ia = {K\f3'\K) = E / (^Va,M^)\"dro , (23) 

KeAi •' 

and the average values of the deformation parameters /5 and 7 in the state \al) are calculated 
from: 

()3);„ = ^fWhc, (24) 

The mixing of different intrinsic configurations in the state \al) can be determined from the 
distribution of the projection K of the angular momentum / on the z axis in the body-fixed 
frame: 

Nk = Q / Ki,(/3,7)l'/?1 sin37M/5c?7, (26) 
Jo Jo 

where the components ipaKH^il) defined in Eq. fl20|) . For large deformations the K 
quantum number is to a good approximation conserved. Consequently, only one of the 
integrals Eq. (1261) will give a value close to 1. A broader distribution of values in the 
state \al) provides a measure of mixing of intrinsic configurations. 

B. Parameters of the collective Hamiltonian 

The entire dynamics of the collective Hamiltonian is governed by the seven functions of 
the intrinsic deformations (3 and 7: the collective potential, the three mass parameters: -B/3/3, 
-B-y-y, and the three moments of inertia X^. These functions are determined by the choice 
of a particular microscopic nuclear energy density functional or effective interaction. As in 
our previous two studies of configuration mixing effects , also in this work we use the 



relativistic functional PC-Fl (point-coupling Lagrangian) [43(] in the particle-hole channel, 
and a density-independent (5-force is the effective interaction in the particle-particle channel. 
The parameters of the PC-Fl functional and the pairing strength constants Vn and Vp have 
been adjusted simultaneously to the nuclear matter equation of state, and to ground-state 
observables (binding energies, charge and diffraction radii, surface thickness and pairing 



gaps) of spherical nuclei [43|], with pairing correlations treated in the BCS approximation. 

The choice of the point-coupling effective Lagrangian determines the self-consistent rela- 
tivistic mean-field energy (RMF) of a nuclear system in terms of local single-nucleon densities 
and currents: 



Ermf = y ^rmf(t') 

= ^ dr vl iJk{r) (-Z7V + m) ^k{r) 
k 

+ -^JTvijTv)f. + —Jtv/^Utv),. + -YPtS + —PTS^PTS + ^PpA 1 , (27) 

where ip denotes the Dirac spinor field of a nucleon. The local isoscalar (S) and isovector 
scalar (TS) densities, and corresponding isoscalar and isovector (TV) currents for a nucleus 
with A nucleons 

psir) = J24Mr)Mr) , (28) 

k 

pTs{r) = ^vl 'ipkir)T3ij^{r) , (29) 

k 

fir) = J2<Mr)YMr) , (30) 

k 

fTvir) = Y.''lMr)rT,^,{r) , (31) 

k 

are calculated in the no-sea approximation: the summation in Eqs. (l27|l - (13T1) runs over 
all occupied states in the Fermi sea, i.e. only occupied single-nucleon states with positive 
energy explicitly contribute to the nucleon densities and currents, vl denotes the occupation 
factors of single-nucleon states. In Eq. ( 1271) pp is the proton density, and denotes the 
Coulomb potential, a, /5, 7 and 5 denote the 11 parameters of the PC-Fl relativistic density 
functional in the corresponding space-isospace channels. 
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The single-nucleon wave functions represent self-consistent solutions of the Dirac equa- 
tion: 

|a ■ HV - V{r)] + V{r) +f3{m + S{r)) } ^,(r) = e,^,(r) . (32) 
The scalar and vector potentials 

S{r) = ^sir) + r^^Ts{r) , (33) 

y^(r) = S^(r) + r3S^^,(r), (34) 

contain the nucleon isoscalar-scalar, isovector-scalar, isoscalar-vector and isovector-vector 
self-energies defined by the following relations: 

= asps + PspI + IspI + Ss^ps , (35) 
St5 = cutsPts + Sts'^Pts , (36) 
= avf + Mjujlf + SvAf - eA^^^ , (37) 

^Tv = (^tvJtv + ^tv^Jtv ' (38) 

respectively. Because of charge conservation, only the 3-rd component of the isovector 
densities and currents contributes to the nucleon self-energies. In this work we only consider 
even-even nuclei, i.e. time-reversal invariance is assumed, which implies that the spatial 
components of the single-nucleon currents vanish in the nuclear ground state. 

The Dirac equation ( l32l) is solved by expanding the nucleon spinors in the basis of a 
three-dimensional harmonic oscillator in Cartesian coordinates. In this way both axial and 
triaxial nuclear shapes can be described. In addition, to reduce the computational task, it 
is assumed that the total densities are symmetric under reflections with respect to all three 
planes xy, xz and yz. When combined with time-reversal invariance, this also implies that 
parity is conserved. Under these restrictions we consider only even-multipole deformations, 
whereas solutions for odd multipoles vanish. The method of solution of the Dirac equation 
is described in more detail in Appendix Rl 

In addition to the self-consistent mean-field potential, for open-shell nuclei pairing cor- 
relations have to be included in the energy functional. In this work pairing is treated using 
the BCS formalism. Following the prescription from Ref. [4^, we employ a 5-force in the 
pairing channel, supplemented with a smooth cut-off determined by a Fermi function in the 
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single-particle energies. The pairing contribution to the total energy is given by 

= / Sl!^hr)dr = .;(„)(r).,(„)(r)rfr , (39) 

for protons and neutrons, respectively. /tp(n)(T') denotes the local part of the pairing tensor, 
and Vp(^n) is the pairing strength parameter. 

The center-of-mass correction is included by adding the expectation value 

to the total energy. Finally, the expression for the total energy reads 

^tot = J [^RMF(r) + S^^M + ^pair(^)] dv + E,^ . (41) 

The entire map of the energy surface as function of the quadrupole deformation is obtained 
by imposing constraints on the axial and triaxial mass quadrupole moments. The method 
of quadratic constraints uses an unrestricted variation of the function 

(^)+ E^2,((Q2,)-g2,)', (42) 

^l=o,2 

where {H) is the total energy, and {Q2fj.) denotes the expectation value of the mass 
quadrupole operator: 

Q20 = - x"" - and Q22 = - y"" . (43) 
q2^ is the constrained value of the multipole moment, and C2fj, the corresponding stiffness 



constant {4]. 

The single-nucleon wave functions, energies and occupation factors, generated from con- 
strained self-consistent solutions of the RMF+BCS equations, provide the microscopic input 
for the parameters of the collective Hamiltonian. 



The moments of inertia are calculated according to the Inglis-Belyaev formula: 2J, |4^ 



x.=5:^^^^l^f^K.i4b)p k= 1,2,3, (44) 

where k denotes the axis of rotation, and the summation runs over the proton and neu- 
tron quasiparticle states. The quasiparticle energies Ei, occupation probabilities Vi, and 
single-nucleon wave functions ipi are determined by solutions of the constrained RMF-I-BCS 
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equations. The mass parameters associated with the two quadrupole collective coordinates 
Qo = {Q20) and q2 = {Q22) are also calculated in the cranking approximation [l8| 



B,,AQo,q2) = Y [■^(i)-^(3)-^(i5 



111/ 



with 



M(n),iii/{q0,q2) 



{i\Q2^\3) {j\Q2u\i) , , x2 



(45) 



(46) 



The collective energy surface includes the energy of the zero-point motion, which has 
to be subtracted. The collective zero-point energy (ZPE) corresponds to a superposition 
of zero-point motion of individual nucleons in the single-nucleon potential. In the general 
case, the ZPE corrections on the potential energy surfaces depend on the deformation. The 
ZPE includes terms originating from the vibrational and rotational kinetic energy, and a 
contribution of potential energy 



^^(90, ^2) = AKib(go, q2) + AVrot(go, ^2) + AVpot(go, q2) ■ 



(47) 



The latter is much smaller than the contribution of kinetic energy, and is usually ne- 



glected 25J. Simple prescriptions for the calculation of vibrational and rotational ZPE 
have been derived in Ref. 18|. Both corrections are calculated in the cranking approxima- 
tion, i.e. on the same level of approximation as the mass parameters and the moments of 
inertia. The vibrational ZPE is given by the expression: 



AKib(go, 92) = ^Tr 



The rotational ZPE is a sum of three terms: 



^(3;a^(2) 



(48) 



AKot(go, ^2) = AV_2-2(go, ^2) + AV_i_i(go, ^2) + AVii(go, ^2), 



(49) 



with 



AV/,^(go,g2) 



1 M(^2),f,u{qo,q2) 



(50) 



4 A<(3),^^(go,g2) ' 

The individual terms are calculated from Eqs. ( l50i) and (H6i) . with the intrinsic compo- 
nents of the quadrupole operator defined by: 



Q21 = -2iyz , Q2-1 = -2xz , Q2-2 = 2ixy 



(51) 
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The potential V^ou in the collective Hamiltonian Eq. (fT2|) is obtained by subtracting the ZPE 
corrections from the total mean-field energy defined in Eq. (j^Tj) : 

Koii(go, ^2) = -E'tot(go, ^2) - AKib(go, 52) - AV^ot(go, ^2) • (52) 

Detailed expressions for the parameters of the collective Hamiltonian are given in Ap- 
pendix [B] and O 

III. ILLUSTRATIVE CALCULATIONS: THE GADOLINIUM ISOTOPIC CHAIN 

In this section the new implementation is tested in a series of illustrative calculations 
of potential energy surfaces and the resulting collective excitation spectra of the chain of 
even-even Gd isotopes: ^^^"^^°Gd. The transition between different shapes, from the weakly- 
deformed transitional ^^^Gd to the well-deformed prolate ^^°Gd, is illustrated in Fig. [H 
where we plot the self-consistent RMF-I-BCS binding energy curves for the axially symmetric 
configurations as functions of the deformation parameter (3. Negative values of (3 correspond 
to the /5 > 0, 7 = 180° axis on the /3 — 7 plane. ^^^Gd is characterized by the coexistence 
of two weakly-deformed prolate {(3 ~ 0.2) and oblate {jS ~ —0.2) minima, with the prolate 
minimum ^ 4 MeV below the spherical configuration. With the addition of more neutrons 
the deformed minima become deeper and gradually shift to larger values of f3. For ^^°Gd, 
the constrained RMF+BCS calculation with the PC-Fl interaction predicts a pronounced 
prolate minimum at {f3 ~ 0.35), more than 10 MeV below the corresponding spherical 
configuration. 

In Figs. [2J12]we display the self-consistent RMF+BCS triaxial quadrupole binding energy 
maps of ^^^"^^''Gd in the /5 — 7 plane (0 < 7 < 60°), obtained by imposing constraints on the 
expectation values of the quadrupole moments (Q20) and {Q22) (cf- Eq. f H2|) ). Filled circles 
denote absolute minima; all energies are normalized with respect to the binding energy of 
the absolute minimum. The contours join points with the same energy. The energy maps 
nicely illustrate the gradual increase of deformation of the prolate minimum with increasing 
number of neutrons. One notes, however, that the oblate minima shown in the axial plots 
in Fig. [1] are not true minima, but rather saddle points in the /3 — 7 plane. 

Starting from constrained self-consistent solutions, i.e. using the single-particle wave 
functions, occupation probabilities, and quasiparticle energies that correspond to each point 
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on the energy surfaces shown in Figs. [2ll6l the parameters that determine the collective 
Hamiltonian: mass parameters -B/3/3, -B/37, -B77, three moments of inertia X^, as well as the 
zero-point energy corrections, are calculated as functions of the deformations (3 and 7. As 
an illustration, for ^^"^Gd the contour map of the inertia parameter (cf. Eq. f lBlQp ) is 
shown in Fig. [TJ and the mass parameter Bpj^ in Fig. [HI The inertia parameter decreases 
with the increase of the axial deformation /3, but the dependence on 7 is very weak. The 
mass parameters, on the other hand, display a pronounced dependence on both intrinsic 
deformations (3 and 7. We notice that Bpp generally increases with the increase of 7 from 
prolate (7 ^ 0°) toward oblate (7 ~ 60*^) configurations. The somewhat erratic behavior of 
Bpp, in particular, is mainly caused by the fluctuations of pairing correlations as function of 
/? and 7 (cf. Eq. (jlHD with its strong dependence on the quasiparticle energies). This effect 
is illustrated in Figs. [9] and [TOl where we plot the contour maps of the proton and neutron 
pairing energies, respectively, in the /? — 7 plane. The fluctuations of pairing energies reflect 
the underlying shell structure and, because pairing correlations are described in the BCS 
approximations, pairing is strongly reduced wherever the level density around the Fermi 
level is small. As a result, mass parameters are locally enhanced in regions of weak pairing. 

Fig. [11] displays the contour plot of the rotational zero point energy correction Eq. fH9|) 
for ^^^Gd. The rotational ZPE, of course, increases with the axial deformation (3, but we also 
notice a pronounced dependence on 7. The ZPE corrections are of the order of several MeV 
even in the region close to the minimum and can, therefore, present a significant contribution 
to the potential of the collective Hamiltonian. For ^^^Gd this is illustrated in Fig. [121 where 
we plot the potential Kou (Eq. (I52l) ) in the /3 — 7 plane. When this plot is compared to 
the total mean- field energy Eq. (14T!) (cf. Fig. [3]), one notes that the main effect of ZPE 
corrections is to shift the position of the minimum to a larger prolate deformation, and to 
modify the shape of the potential around the minimum. 

The diagonalization of the resulting Hamiltonian yields the excitation energies and the 
collective wave functions for each value of the total angular momentum and parity I'" . 
In addition to the yrast ground-state band, in deformed and transitional nuclei excited 
states are usually also assigned to (quasi) (3 and 7 bands. This is done according to the 
distribution of the angular momentum projection K quantum number defined in Eq. (l26l) . 
Excited states with predominant K = 2 components in the wave function are assigned to 
the 7-band, whereas the /3-band comprises states above the yrast characterized by dominant 
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K = components. As an example, in Fig. [T3l w e display the PC-Fl excitation spectrum 
of ^^^Gd, in comparison with available data [45]. The relative excitation energies within 
all three theoretical bands are scaled by the common factor ^ 0.69, determined in such a 
way that the calculated energy of the 2f state coincides with the experimental value. This 
leaves the bandheads of the 7 and (3 bands unaltered, and corresponds to an enhancement of 
the effective moment of inertia by ~ 45%. The scaling of the relative excitation energies is 
introduced because of the well known fact that the Inglis-Belyaev (IB) formula (jH]) predicts 
effective moments of inertia that are considerably smaller than empirical values. More 
realistic values are only obtained if one uses the Thouless-Valatin (TV) formula {22], but 
this procedure is computationally much more demanding, and it has not been implemented 
in the current version of the model. Here we rather follow the prescription of Ref. |25j] where, 
by comparing the TV and IB moments of inertia as functions of the axial deformation for 
superdeformed bands in the A = 190 — 198 mass region, it was shown that the Thouless- 
Valatin correction to the perturbative expression IB is almost independent of deformation, 
and does not include significant new structures in the moments of inertia. It was thus 
suggested that the moments of inertia to be used in the collective Hamiltonian can be simply 
related to the IB values through the minimal prescription: Xjt(g) = Tl^{q)(l + a), where q 
denotes the generic deformation parameter, and a is a constant that can be determined in 
a comparison with data. The value of a ~ 0.45 used for the excitation spectrum of ^^^Gd is 
somewhat larger than the values determined in the mass A = 190 — 198 region 251]. 

When the IB effective moments of inertia are renormalized to the empirical values by 
scaling the relative excitation energies to reproduce the experimental position of the state 
2'1, the resulting bands are in good agreement with data. This is illustrated in Figs. [T^lfT6l 
where we plot the excitation energies with respect to bandheads, for the ground-state band, 
7, and (3 bands in ^^^^^^'^Gd, respectively. For each nucleus the relative excitation energies 
within the three bands are scaled by a common factor, adjusted to the experimental energy 
of the 2f state, as explained above. These factors are rather similar for four isotopes: 0.69 
for is^Gd, 0.67 for ^^^Gd, 0.69 for ^^^Gd, and 0.72 for ^^^Gd. An exception is the lightest 
Gd isotope considered here: ^^^Gd, for which this factor is actually 1.08, i.e. the theoretical 
2f state {E^ = 0.318 MeV) is slightly below the experimental 2f state {E^ = 0.344 MeV). 
However, as shown in Figs. [T4lfT6l for this nucleus the calculated ground-state band, as 
well as the (quasi) 7 and /3 bands, do not follow very closely the experimental spectra. The 
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deviation from the experimental trend at higher angular momenta can probably be explained 
by the fact that this is a weakly-deformed transitional nucleus for which the assumption of 
a constant moment of inertia, implicit in the expression Xfc(g) = X^^(g)(l + a), and of 
relatively pure (3 and 7 bands, does not present a very good approximation. 

The calculated (3 bands are compared with available data in Fig. [TSl In comparison with 
the 7 bands (cf. Fig. [T5|l . the agreement with data is better in ^^^Gd, but the deviation from 
experiment is more pronounced in ^^^Gd. To understand in more detail the discrepancy 
between the calculated and empirical f3 and 7 bands, we need to consider the distribution 
of the angular momentum projection K quantum number (cf. Eq. ( 126|) ) in these bands. 
As explained above, the calculated second and third eigenstate for each angular momentum 
are assigned either to the /3 or 7 band, on the basis of the predominant K = or K = 2 
components, respectively. The distributions of K components in the wave functions of the 
calculated second and third 2+, 4+ and 6'*' states are plotted in Figs. [T71fT9| respectively. 
In the case of ^^^Gd, in particular, we find a pronounced mixing of the K = or K = 2 
components and, with increasing angular momentum, contributions of higher- 7^ components 
are present in the wave functions. This is consistent with the observation that ^^^Gd is a 
transitional nucleus and, therefore, excited states can only approximately be assigned to 
(quasi) f3 and 7 bands. Increasing the neutron number toward heavier and more deformed 
Gd isotopes, the distributions of K components become sharp and states can unambiguously 
be grouped into (3 and 7 bands. This is, of course, characteristic for well deformed nuclei. 
One exception, however, is the calculated spectrum of ^^^Gd, where we find significant mixing 
of i^' = and K = 2 components in the wave functions of second and third 2"*", 4"'" and 6"*" 
states, as well as for higher angular momenta. The more pronounced mixing between the 
P and 7 bands occurs because, in this particular isotope, the calculated second and third 
even-spin states are almost degenerate in energy, as shown in Fig. [201 

The level of i^'-mixing is reflected in the staggering in energy between odd- and even-spin 
states in the (quasi) 7-bands (cf. Fig. [15]). The staggering can be quantified by considering 
the differential quantity 



{i?[J+] - E[iJ - 1)+]} -{E[{J - 1)+] - E[iJ - 2)+]} 

E[2t 



g^j-^ ^ X^L^^j ^LV^ x^iK'^ ^jyi ^LV^ ^^^^ 



which measures the displacement of the [J — l)^ level relative to the average of its neighbors, 
and (J — 2)+, normalized to the energy of the first excited state of the ground state 
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band, 2^. Because of its differential form, S{J) is very sensitive to structural changes. For 
an axially symmetric rotor S{J) is, of course, constant. In a nucleus with a deformed 7-soft 
potential, S{J) oscillates between negative values for even-spin states and positive values for 
odd-spin states, with the magnitude slowly increasing with spin. For a triaxial potential the 
level clustering in the (quasi) 7-band is opposite, and S{J) oscillates between positive values 
for even-spin states and negative values for odd-spin states. In this case the magnitude of 
S{J) increases more rapidly with spin, as compared to the 7-soft potential. In a recent study 
of staggering of 7-band energies and the transition between different structural symmetries 



in nuclei 



471], the experimental energy staggering in 7-bands of several isotopic chains has 



been investigated as a signature for the 7 dependence of the potential. In Fig. HT] we plot 
the staggering in the 7-band for the chain of the Gd isotopes, calculated with the PC-Fl 
relativistic density functional. One notices how the pronounced i^-mixing in ^^^Gd and 
^^^Gd (cf. Figs. [T71fT9|) leads to the strong staggering observed in the corresponding (quasi) 
7-bands. The calculation reproduces both the empirical oscillatory behavior and, with the 
exception of low-spin states in ^^^Gd, also the magnitude of S{J). Starting from the 7-soft 
^^^Gd (negative values for even-spin states and positive values for odd-spin states), S{J) 
evolves toward the axially symmetric rotor limit {S{J) = 0.33) in ^^^Gd and ^^°Gd. 

The assignment of even-spin states above the yrast either to the /3 or 7 bands, on the basis 
of the predominant K = or K = 2 components (Figs. [T71fT9|) . and the level of K- mixing 
inferred from the differential quantity S{J) Eq. (|53l) (Fig. [2T1) . has a correspondence in the 
calculated average values of the deformation parameters f3 and 7 (cf. Eqs. ( 12^ and f l25l) ). 
In Tabled] we collect the average /3 and 7 deformations for the calculated first, second, and 
third 2"*", 4"*" and 6"^ states in ^^^~^^°Gd. For those nuclei where the K-mixing is weak (sharp 
distribution of K-components in Figs. [T71fT9| weak staggering of S{J) in Fig. [21]): ^^''Gd, 
^^^Gd, and ^^'^Gd, the average values of the deformation parameters are almost identical for 
states belonging to the ground-state band and those assigned to the /?-band. States assigned 
to the 7-band {K = 2) are consistently characterized by larger average values of the angle 
7. This distinction does not appear in the spectra of the two nuclei for which the model 
predicts pronounced K-mixing: ^^^Gd and ^^^Gd. 

An important advantage of using structure models based on self-consistent mean-field 
single-particle solutions is the fact that physical observables, such as transition probabilities 
and spectroscopic quadrupole moments, are calculated in the full configuration space and 
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there is no need for effective charges. Using the bare value of the proton charge in the 
electric quadrupole operator j(4{E2), the transition probabilities between eigenvectors of 
the collective Hamiltonian can be directly compared with data. In addition to the calculated 
energy spectrum of ^^^Gd, in Fig. [12] we have also compared the resulting B(E2) values (in 
Weisskopf units) for transitions within the ground-state band, and the transitions 2^ — > 0^ 
and 0^ — i> 2g g , with available experimental values. The agreement between theoretical 
B(E2) values and data is very good, especially considering that the calculation of transition 
probabilities is completely parameter-free. We also notice the remarkable prediction for 
the interband transition 0^ — > 2^^ , in excellent agreement with experiment. Finally, in 
Fig. [22] we plot the calculated B{E2) values (in Weisskopf units) for the ground-state band 
transitions ^ {J ~ 2)^ in ^^^"^^'^Gd, together with the available experimental values. 
The model clearly reproduces the empirical trend of ground-state band transitions in Gd 
isotopes and, except perhaps for the transitional nucleus ^^^Gd, the theoretical predictions 
are in excellent agreement with data even for higher angular momentum states. 



IV. SUMMARY AND OUTLOOK 



To describe complex excitation patterns and electromagnetic transition rates associated 
with the evolution of shell structures starting from stable nuclei and extending toward re- 
gions of exotic short-lived systems far from /5-stability, nuclear structure methods must be 
developed that are based on a universal microscopic framework. Properties of a vast major- 
ity of nuclides with a large number of valence nucleons are best described by nuclear energy 
density functionals. However, for a quantitative description of energy spectra and transition 
probabilities one must be able to go beyond the lowest order in which the EDFs are imple- 
mented - the mean-field approximation, and systematically include correlations related to 
the restoration of broken symmetries and to fluctuations of collective variables. In the frame- 
work of non-relativistic EDFs, in recent years several models have been developed that use 
the generator coordinate method (GCM) to perform configuration mixing calculations with 
angular-momentum and particle-number projected mean-field (SR EDF) states. In most 
applications the calculations have been restricted to axially symmetric, parity conserving 
configurations. 

In a recent series of papers, of which the present is the third part, we have extended 
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the relativistic EDFs to include symmetry restoration and fluctuations of quadrupole de- 
formations. While in the first two parts the GCM was used in configuration mixing cal- 
culations with axially symmetric relativistic wave functions, this work has been focused 
on the description of general triaxial shapes. We have developed a new implementation 
for the solution of the eigenvalue problem of a five-dimensional collective Hamiltonian for 
quadrupole vibrational and rotational degrees of freedom, with parameters determined by 
constrained self-consistent relativistic mean-field calculations for triaxial shapes. In addition 
to the self-consistent mean-field potential of the PC-Fl relativistic density functional in the 
particle-hole channel, for open-shell nuclei pairing correlations are included in the BCS ap- 
proximation. The resulting single-nucleon wave functions, energies and occupation factors, 
as functions of the quadrupole deformations, provide the microscopic input for the param- 
eters of the collective Hamiltonian: three mass parameters: -B/3/3, -B/37, -B77, three moments 
of inertia X^, and the collective potential including zero-point vibrational and rotational 
energy corrections. The moments of inertia are calculated using the Inglis-Belyaev formula, 
and the mass parameters associated with the quadrupole collective coordinates are deter- 
mined in the cranking approximation. An extensive test has been carried out in calculations 
of potential energy surfaces, and the resulting collective excitation spectra and transition 
probabilities, for the chain of even-even gadolinium isotopes. Results for excitation energies 
in ground-state, (quasi) (3 and 7 bands, and the corresponding interband and intraband 
transition probabilities have been compared with available data on even-even Gd isotopes: 

152-160Qd_ 

There are some obvious improvements that need to be implemented in the model. For 
instance, because the Inglis-Belyaev formula gives effective moments of inertia that are lower 
than empirical values, all the calculated relative excitation energies in ^^^"^^''Gd had to be 
scaled with respect to the experimental energy of the 2f states. The moments of inertia can 
be improved by including the Thouless-Valatin dynamical rearrangement contributions. For 
the rotational degrees of freedom for which the collective momenta are known, the inertia 
parameters can be obtained from the solutions of cranked RMF equations. For the defor- 
mation coordinates go and q2 the situation is more complicated, because the corresponding 
momentum operators Pq and P2 have to be calculated from the solution of Thouless-Valatin 



equations 



22j at each deformation point. Because cranking breaks time- reversal symmetry. 



in both cases the inclusion of pairing correlations necessitates the extension of the model 
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from the simple RMF+BCS to the full relativistic Hartree-Bogoliubov framework 



APPENDIX A: THREE-DIMENSIONAL SOLUTION OF THE DIRAC EQUA- 
TION 

To solve the Dirac equation (!32l) for triaxially deformed potentials, the single-nucleon 



spinors are expanded in the basis o: 
(HO) in Cartesian coordinates 48 



eigenfunctions of a three-dimensional harmonic oscillator 
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50, l5l(| • In one dimension: 



= b;'/'H^^{Qe-^'^/' , (/i ^ a;, y, z) (Al) 
= Xfj_/b^, and the oscillator length is defined as 



K = xl-^- (A2) 



HniO denotes the normalized Hermite polynomials 

/oo 
-oo 



(A3) 



The basis state can be defined as the product of three HO wave functions (one for each 
dimension) and the spinor: 

$a(r;m,) = 0n,(^x)0n,(^s,)0n,(^^)Xm,, (A4) 

where the notation is: a = {ux, ny,nz}. It will be more convenient to use the eigenstates of 
the x-simplex operator defined by the relation 

Sx = Pe-'''^% (A5) 

where P denotes the parity operator. It is easily verified that the x-simplex operator acting 
on the state $a(r;ms) leads to 

4$a(r; m,) = -z(-l)"-<l>„(a;, y, z; -m,) . (A6) 

The eigenstates of the x-simplex operator with positive and negative eigenvalues read: 

$„(r; +) = 0„^(a;)0„^(y)0„^(z)^ [x+ - (-l)"-^-] (A7) 

<|.„(r; -) = 0„^(x)0„^(y)0„,(^)^(-l)-+"«+i [X+ + i-irx^] . (A8) 
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For the Dirac spinor with positive simplex eigenvalue, the large component corresponds 
to positive, and the small component to negative eigenvalues 




^^{r,+)= (A9) 
The large and small component are expanded in terms of the basis states Eqs. ( IA7I) and 

h{v- +) = 5^ /r$„(r; +) and g,{v- -) = ") " (^1°) 

a a 

To avoid the occurrence of spurious states, amax and oimax are chosen in such a way that the 
corresponding major quantum numbers N = nx + Uy + Uz are not larger than some arbitrary 
Np for the expansion of large components, and not larger than Np + 1 for the expansion of 
small components ^]. The single-nucleon Dirac equation: 



V + m* (T-V \ fi{r; 
~<7 ■ V V -m* I \ gi{r; 




(All) 



with the effective nucleon mass m* = 771 + 3, reduces to a symmetric matrix eigenvalue 
problem 

Aaa' Baa' \ | A" ) _ | fi 
Baa' Caa' ) \ 9?' ) \ 9? 



(A12) 



of dimension: amax + ctmax- 

The time-reversal operator T = —lOyK exchanges the simplex eigenvalues: 

T<l>„(r;+) = -<l>„(r;-) and f $„(r; -) = <l>„(r; +) , (A13) 

and, thus, when acting on the Dirac spinors 

T^.(r;+)=T| 1 = 1 -r^' 1 = 1 1 = -^..(r; -) . (AM) 



/ /.(r;+) 1 


H 


V Wiij; -) j 






-i9iir; +) 

Time-reversed single particle states correspond to opposite simplex eigenvalues. Because of 
time-reversal invariance, for each solution of the Dirac equation flAlip with positive simplex 
eigenvalue 7pi{r; +), there exists a degenerate time-reversed solution with negative simplex 
eigenvalue ipi^r; —). Both solutions contribute equally to the densities, and in practice only 
the Dirac equation for positive simplex eigenstates is solved. 
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In the current implementation of the model, parity is also imposed as a conserved symme- 
try. This means that the basis states split into two parity blocks, which can be diagonalized 
separately. In addition, we require that the densities are symmetric under reflections with 
respect to the yz^ xz and xy planes. 

ps,v{.x, y, z) = ps,vi-x, y, z) = ps,vix, -y, z) = ps,vix, y, -z) . (A15) 

The symmetries of the scalar and vector densities are, of course, fulflUed by the corresponding 
self-consistent scalar and vector potentials: 

S{,x, y, z) = S{-x, y, z) = S{x, -y, z) = S{x, y, -z) (A16) 
Vix, y, z) = V{-x, y, z) = V{x, ~y, z) = V{x, y, -z) . (A17) 

The self-consistent symmetries (]A16P and (]A17|) simplify the evaluation of the matrix ele- 
ments Aaa' and Caa'- First, the symmetry under reflections with respect to the xy, yz and 
xz planes means that we need to calculate only matrix elements between states (pa and (pa' 
for which: 

(-1)"- = (-l)< , (-1)"^ = (-1)"; and (-1)"^ = (-1)'^^ . (A18) 

Furthermore, three-dimensional integrals reduce to the octant x,y,z > 0. The matrix 
elements of the vector and scalar potentials read: 

^aa' = + r (A19) 

POO poo poo 

Jo Jo Jo 

Caa' = {a;~\V -m*\a';-) (A20) 

POO poo poo 

= S{-lY-y-<y^ / / / 0„Jx)0„^(i/)0„^(^)(l^-m*)0„,(x)0„;(y)0„,(^)c?V^. 
Jo Jo Jo 

Note that the condition (lAlSP means that the difference [uy — n' ) must be even, hence the 



matrix elements (1A19P and (]A20p are real. 



The matrix elements of the kinetic energy term Baa' can be calculated analytically using 
the expression: 

dt,K^{xt,) = [-^/n~rT(l)n^+i{x^) + v^0n^-i(a;/.)] , (A21) 
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together with the orthogonahty relation (lA3p . The operator a ■ V consists of three terms: 



Ba&' = (a; + \-a^d^ - aydy - a^d^\ a'; -) = B^^, + Bl^, + B^^, 



(A22) 



and the corresponding matrix elements are calculated from 

1 



Baa' — ( '^T^ ^nyn'y^n^rv^ 
BL, = {-ir'^5n.K5n.n' 



_ /I N<+n' +1 r r 



V2h 



n'y + + JK^nyn'y-l 



(A23) 
(A24) 
(A25) 



The set of self-consistent solutions of the single-nucleon Dirac equation determines the 
scalar and vector densities: 



p,(r) = 2 5^t;2^1(r;+)^,(r;+), 

i 



(A26) 
(A27) 



Because of time-reversal symmetry, the summation is over positive simplex solutions. To 
calculate densities in coordinate space, one needs explicit expressions for the products of 
basis states: 

0t (r; +)0,,(r; +) = (-l)K-«.)/20„^(a;)0„, (^), (A28) 

0l(r; -)0«Kr; -) = (-l)("^-"^)/'</)„.(x)0„, (^). (A29) 



Note that the symmetry requirement Eq. (lAlSp imposes the condition (lAlSp and, again, the 
difference (n^ — n'y) is even so that the contributions (]A28|) and (]A29|) to the densities are 
both real. 



APPENDIX B: MOMENTS OF INERTIA 

The basic ingredient of the Inglis-Belyaev formula for the moments of inertia Eq. (jH]) 
are the matrix elements of the angular momentum operators in the simplex basis (]A7[1 and 
( lASp . Here we present in detail the calculation of the matrix element of the Jx component 
between basis states with positive simplex eigenvalue 

. (Bl) 
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For the other matrix elements only the final expressions will be listed. 

The X component of the total angular momentum operator is the sum of the spin and 
the spatial contributions: 

h - h 

Jx = -^^x + L^ = -a^ - ih {ydz - zdy) . (B2) 

The spatial parts of the basis states are unaffected by the dx operator, thus generating the 
product of Kronecker delta symbols 5n^^n'^5ny^n'y5nz,n'-^- The contribution from the spin factors 
of the basis states is given by 

[xl - - {-IT-XA = i-ir^' + (-1)"^+' , (B3) 

and the spin matrix element reads: 

^{a; +\ax\a'; +) = ^5„^,<5„„„;,5„,,„i(-l)"^+' . (B4) 
Next, the contribution from the operator is calculated 

(a; +\Lx\a'; +) = -ih{a; +\ydz - zdy\a'; +). (B5) 
The spin factors of the basis states are not affected by the operator: 

[xl - - i-l^X-] = 1 + (B6) 

To calculate the matrix elements Eq. (IBSp . the following relations are used: 

a^M^n^ = ^ [^/ni^TT(|)n^+l{x^) + y^(j)n^-i{x^)] , (B7) 



2 "-"^ \h by 

by 



d^(l)n^iXf,) = — =— [-^/Uf, + + yrv0„,^-i(a;^)] , (B8) 

together with the orthonormality relation (]A3|) . The total matrix element reads 
(a; +\l\a'; +) = 5„^_„, (B9) 

+ 7:^n^.< (t^ - 7^1 ^ra^ + 1a/< + 16ny,n'y+lSn,,n',+l + \J^y^/^z^ny,n'y-l8n,,n', 

Uy ,n' —l^Uz ,n'^ 

The following relations can easily be proved: 

(a; +1 J^^la'; +) = — (a; — I Ja^la'; — ), and (a; +| Ja;.|a'; — ) = (a; — | JjIq;'; +) = . (BIO) 
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+1 



The final expression for the moment of inertia Xr^ =Ii: 

{UiVj - ViUjf 



i,j>0 



Ei + Ej 



(Bll) 



We include only the final results for the matrix elements of Jy and J^: 

h 

{a; +\Jy\a'; -) = ^-(-l)"''5n.,ni^"j/-<^"..< 



x "'y>"'y 



(B12) 



(a; -| Jj^la'; +) = -(a; +1 Jj^la'; -), and (a; +| Jj^|q;'; +) = (a; -| Jj^|q;'; -) = 0, (B13) 



(B14) 



6„ 6^ 



X 



(a;-|J,|a';+) = (a;+|J,|a';-), and (a; +| J,|a'; +) = (a; -| J,|a'; -) = 0, (B15) 
The corresponding moments of inertia Xy = X2 and = X^ read: 



i,j>0 



-Ej- + Ea 



aa' aa' 



, (BI6) 



(BIT) 



All three moments of inertia vanish at the spherical point (3 — 0. In addition, X^ vanishes 
for the 7 = 0° configurations (prolate deformed, with z as the symmetry axis), whereas 
Xy vanishes at 7 = 60*^ (oblate deformed, y is the symmetry axis). These conditions are 
incorporated in the following functional form: 



Xk = 4Sfe/3^ sin^ (7 - 2kn/3) , with (1 = x, 2 = y, 3 = z) , 



(B18) 
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from which the inertia parameters B^, By and Bz follow: 



Bk 



2fc 



(B19) 



sin^ (7 - 2kn/3)' 



For the limiting cases described above the following relations are used 



B,{P = 0) = ByiP = 0) = Bz{/3 = 0) = B^^iP = 0) 



(B20) 



5, (A 7 = 60°) =5^^ (A 7 = 60°), 



(B21) 



5,(^7 = 0°) = 5,,(A 7 = 0°). 



(B22) 



To calculate the mass parameters from Eqs. (145|) and (l46l) . one needs the matrix elements 
of the operators x^, y"^ and 2;^ in the simplex basis. These are combined and inserted into 
the matrix elements 



m Q2, 1^,) = Yl ftff{<^-, +\Q2,\a'; +) + J29f9f{^; +\Q2,\a; +) , (B23) 



that determine the 2x2 matrix M.{>n),^v Eq. P6|) . The mass parameters B^^{qo,q2) are 
then calculated from Eq. (HSl) . and transformed from the quadrupole coordinates qo,q2 to 
the polar Bohr deformation variables (3 and 7. 

APPENDIX C: ROTATIONAL ZERO-POINT ENERGY CORRECTION 

The rotational zero-point energy Eq. (149|) is determined by the matrix elements of the 
quadrupole operators: 



Q21 = -2iyz , Q2-1 



2xz , and Q2-2 = 2ixj/ . 



(CI) 



Using the following expression: 




(C2) 
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together with the orthogonahty relation (lA3p . the calculation of matrix elements is straight- 
forward. Here we list only the final expressions: 



X 



X 



(a; +|(52i|a'; +) = bybJn^^n'^ 

-(a; +|Q2i|tt'; +) 
{a;-\Q2i\a';+) = 



{a; -\Q2i\<y'; - 
{a; +\Q2i\a'; - 
(a; +\Q2^i\<y'; - 



{a; -\Q2-1\a'; + 
{a; +\Q2-i\a'; + 
{a; +\Q2-2\a'; - 

{a; -\Q2-2\a'; - 
{a; +\Q2-2\a'; - 



n'-l 



X 



X 



-{a; +\Q2-i\a'; -) 
{a; -\Q2-1\a'; -) = 



n'y + lSny,n'y+l - 

-{a; +\Q2~2\oi'; +) 
{a; -\Q2-2\a'; +) = 



nySny,n'y-l 



(C3) 
(C4) 
(C5) 

(C6) 
(C7) 
(C8) 

(C9) 
(CIO) 
(Cll) 



APPENDIX D: NUMERICAL DETAILS 

In Fig. [23] we check the convergence of the RMF+BCS quadrupole constrained calcula- 
tions as a function of the maximal number of oscillator shells used in the expansion of the 
Dirac spinors. The binding energy curves for ^^^Gd are plotted as functions of the axial (3 
deformation. These curves correspond to calculations with 10, 12, 14 and 16 major oscillator 
shells. For moderate deformations considered in this study {\(3\ < 0.65), the RMF results 
show convergence at Nf = 14. 

The self-consistent RMF+BCS equations are solved on a mesh over the first sextant of 
the /3 — 7 plane: 

/? > 0; Ap = 0.05; < 7 < 60°; A7 = 6° , (Dl) 

and the maximum value of the (3 deformation is: jSmax = 1.15. 

The choice of the basis wave functions Eq. ( fT6|) introduces a factor e~^^^^ into the integral 
over (3 in the matrix elements of the collective Hamiltonian. With the substitution y = 
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we obtain the weight function appropriate for Gauss-Laguerre quadrature. The integrals over 
7 are evaluated by Gauss- Legendre quadrature. The corresponding number of mesh points 
are: = 24 and = 24, respectively. The parameters of the collective Hamiltonian 
in the Gaussian mesh-points are calculated by interpolation from the values calculated on 
the equidistant mesh Eq. (IDip . To avoid extrapolation, the minimum value of the basis 
parameter fi is restricted to: 



We have verified that both the calculated excitation spectra and transition probabilities 
remain stable for any choice of fi within the interval: 8 < yU < 15. All calculations presented 
in this work are performed with /i = 9. 
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FIG. 1: (Color online) Self-consistent RMF+BCS binding energy curves of the ^^^"^^'^Gd isotopes, 
as functions of the axial deformation parameter /3. Negative values of /3 correspond to the {f3 > 
0, 7 = 180*^) axis on the /3 — 7 plane. 
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FIG. 2: (Color online) Self-consistent RMF+BCS triaxial quadrupole binding energy map 0fl52Gd 
in the /? — 7 plane (0 < 7 < 6O''). All energies are normalized with respect to the binding energy 
of the absolute minimum (red dot). The contours join points on the surface with the same energy 
(in MeV). 
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FIG. 3: (Color online) Same as Fig. [21 but for the nucleus ^^^Gd. 
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FIG. 4: (Color online) Same as Fig. [21 but for the nucleus ^^^Gd. 
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FIG. 5: (Color online) Same as Fig. [21 but for the nucleus ^^^Gd. 
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FIG. 6: (Color online) Same as Fig. [21 but for the nucleus ^^''Gd. 
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FIG. 7: (Color online) The map of the Inglis-Belyaev inertia parameter Bx of ^^''Gd in the /? — 7 
plane (0 < 7 < 60°). 
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FIG. 9: (Color online) Proton pairing energy of i^^Gd in the /? - 7 plane (0 < 7 < 60°). The 
contours join points on the surface with the same energy (in MeV). 
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FIG. 10: (Color online) Neutron pairing energy of ^^^Gd in the /3 — 7 plane (0 < 7 < 60'^). The 
contours join points on the surface with the same energy (in MeV). 
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FIG. 11: (Color online) The rotational zero-point energy of ^^^Gd in the /? — 7 plane (0 < 7 < 60*^). 
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FIG. 12: (Color online) The potential VcoW (Eq. ^) of ^^^Gd in the /? - 7 plane (0 < 7 < 60°). 
The contours join points on the surface with the same energy (in MeV). 
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FIG. 13: (Color online) The level scheme of Gd calculated with the PC-Fl relativistic density 



functional, in comparison with the experimental data [45(. The relative excitation energies are 
scaled by the common factor ~ 0.69, adjusted to the experimental energy of the state 2^. The 
B{E2) values are given in Weisskopf units. 
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FIG. 14: (Color online) Relative ground-state band excitation energies in ^^^^^^''Gd. For each 
nucleus the theoretical energies are scaled by a common factor, adjusted to the experimental 
energy of the 2f state. 
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FIG. 17: (Color online) Distribution of the X-components (projection of the angular momentum 
on the body-fixed symmetry axis) in the collective wave functions for the states: 2^ and 2^. 
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FIG. 18: (Color online) Same as in Fig. Wl\ but for the states: 4^ and 43 
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FIG. 20: (Color online) Excitation energies of the second and third states: 2^, 4^ and 6^ in the 
Gd isotopic chain, as functions of the neutron number. 
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FIG. 22: (Color online) B{E2) values (in Weisskopf units) for the ground-state band transitions 
Ji ^ {J — 2)^ in ^^^~^^'^Gd. Theoretical values calculated with the PC-Fl relativistic density 
functional are compared with data. 
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FIG. 23: (Color online) Self-consistent RMF+BCS binding energy curves of ^^^Gd, as functions of 
the axial deformation parameter /3. Negative values of /? correspond to the (/3 > 0, 7 = 180*^) axis 
on the /3 — 7 plane. The four energy curves correspond to calculations in the three-dimensional 
harmonic oscillator basis with 10, 12, 14 and 16 major oscillator shells, respectively. 
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TABLE I: Average values of the deformation parameters f3 and 7 (cf. Eqs. (j24p and ()25p ) for the 
calculated first, second, and third 2+, 4+ and 6+ states in -'^^^"-'^^''Gd. 



2+ 4+ 6+ 





state 




(7) (deg) 




(7) (deg) 




(7) (deg) 






0.24 


17.0 


0.26 


15.3 


0.28 


13.7 


152Gd 


4 


0.26 


19.0 


0.31 


13.9 


0.32 


12.8 




4 


0.29 


16.3 


0.28 


20.5 


0.30 


20.1 




Jt 


0.31 


12.8 


0.32 


12.0 


0.33 


11.3 


154Gd 
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0.33 


12.8 


0.33 


12.4 


0.34 


12.2 
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0.29 


19.8 


0.32 


17.6 


0.34 


16.2 
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0.34 


11.3 


0.35 


11.0 


0.36 


10.6 


156Gd 
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0.34 
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0.36 
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0.36 
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0.36 
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13.4 




4 


0.36 


11.0 


0.37 
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0.38 


10.4 




4 


0.36 


10.3 


0.37 


10.2 


0.37 


10.1 


160 Gd 


4 


0.37 


13.4 


0.37 


13.2 


0.37 


12.9 
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0.38 
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0.38 
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0.39 


9.8 
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